Method for Estimating the Heat Load Imposed on a Cryogenic Refrigerator, Associated Program Product, and Method for Controlling the Refrigerator

ABSTRACT

The invention relates to a method for estimating a heat load imposed on a cryogenic refrigerator, to an associated computer program product, and to a method for controlling the cooling power output by said refrigerator. As the refrigerator ( 1, 1 ′) includes a phase separator ( 40, 40 ′) comprising a bath ( 41, 41 ′) of refrigerant, the method for estimating the heat load imposed on said refrigerator includes a step in which said heat load is estimated using a program executed by a computer, said program being based on a mass balance carried out on the phase separator for expressing variations in the time drift of the height of the bath of refrigerant in the phase separator.

The present invention relates to the field of plant refrigeration.

The present invention in particular relates to the cryogenicrefrigeration of plants.

The present invention also relates to the cryogenic refrigeration ofplants capable of operating in a variable regime.

A variable operating regime may be encountered in many applications.

This is for example the case in plants using strong magnetic fields.

An example of a plant using strong magnetic fields is a tokamak. Atokamak is a magnetic confinement chamber intended to control a plasmain order to study the possibility of power generation by nuclear fusion.

In tokamaks, a variable, pulsed operating regime may be employed. Inthis case, although the overall profile of the pulsed loads applied tothe refrigerator is known, it is not exactly known when a pulse of loadwill occur. In addition, other unforeseeable perturbations may occur,which perturbations are related to the operation of the tokamak.

To generate strong magnetic fields without destroying theelectromagnets, superconducting electromagnets are used. For anelectromagnet to operate as a superconductor, its temperature must bekept below its critical temperature.

This is achieved by virtue of a cryogenic refrigerator, for examplecooled by helium.

FIG. 1 is a schematic showing the operation of such a cryogenicrefrigerator.

The cryogenic refrigerator 1 comprises a compressor 10 allowing a gas,in this case helium, at room temperature (T₀≈300 K) to be compressedfrom atmospheric pressure P_(a) to a pressure P of about 20 bar.

A number of heat exchangers 20, 21, 22, 23, 24 are placed in parallelupstream of a phase separator 40 comprising a bath of liquid helium at atemperature T_(f)=4.5 K. In this case, five heat exchangers have beenprovided.

These countercurrent heat exchangers allow the temperature of the heliumflowing in the various circuits of these heat exchangers to be graduallydecreased. Moreover, for each heat exchanger, the pressure in the twoexchange circuits of the exchanger differ, so that the “upstream”circuit is a hot, high-pressure circuit and the “downstream” circuit isa colder, low-pressure circuit. Between the “upstream” circuit of thefirst heat exchange 20 and the “downstream” circuit of the last heatexchanger 24, the pressure thus increases from P=20 bar at the outlet ofthe compressor to a pressure slightly higher than atmospheric pressure,to prevent problems with low-pressure cavitation.

The bath 41 of liquid helium is therefore at atmospheric pressure.

The cryogenic refrigerator 1 also comprises a number of means forextracting work.

In the present case, these means consist of two turbines 30, 31.

The first turbine 30 has work done on it at the outlet of thelow-pressure circuit of the first heat exchanger 20 and reinjects thiswork at the low-pressure inlet of the second heat exchanger 21. Thesecond turbine 31 has work done on it at the outlet of the low-pressurecircuit of the third heat exchanger 22 and reinjects this work at thelow-pressure inlet of the fourth heat exchanger 23.

These turbines 30, 31 are complementary to the heat exchangers andparticipate, via the work done on them, to the cooling of the helium.

Lastly, the cryogenic refrigerator 1 comprises a Joule-Thomson valve 50,placed between the outlet of the low-pressure circuit of the last heatexchanger 24 and the bath 41 of liquid helium at 4.5 K and atmosphericpressure. This valve 50 liquefies the gaseous helium obtained at theoutlet of the low-pressure circuit of the last heat exchanger 24 via anexpansion which is accompanied by a drop in the temperature of thehelium.

The bath 41 of liquid helium then allows cooling power to be deliveredin order to keep the electromagnets of the plant employing strongmagnetic fields, for example a tokamak, operating as superconductors, onaccount of the heating power applied by the plant to this bath 41 ofliquid helium.

The heating power applied by the plant to the bath 41 of liquid heliumis also called the thermal load.

The design and dimensions of this type of refrigerator are tailored tooperating regimes in which the strong magnetic field generated by theplant is stable or varies slowly, i.e. to a regime of permanent oralmost permanent thermal operation of the refrigerator.

Specifically, these plant operating regimes lead to a stable thermalload on the cryogenic refrigerator. Operating the cryogenic refrigeratorin a permanent or almost permanent thermal regime makes it possible toprovide adequate cooling of the electromagnets, stably and reliably.

However, future plants intended for studying the possibility ofgenerating power by nuclear fusion activated by strong magnetic fieldsplan to employ pulsed magnetic fields.

This is the case for the ITER (France) and JT60_SA (Japan) tokamakprojects.

Variations in the magnetic field in the tokamak then cause similarvariations in the thermal load applied to the cryogenic refrigerator.

An example of the variation in the thermal load applied to the cryogenicrefrigerator of the future Japanese tokamak JT60_SA, intended to operatein a pulsed regime, is shown in FIG. 2. It will be noted that the pulsesare a priori periodic, in this instance having a period of 3000 s, theduration of a pulse being about 100 s. Nevertheless, random processesmay change this variation, as was mentioned above.

In FIG. 2, it may be seen that operating the tokamak in a pulsed regimeleads to the operating regime of the thermal load applied to thecryogenic refrigerator also being pulsed (curve 100). FIG. 2 also shows(curve 101) the average thermal load applied by the tokamak to therefrigerator.

However, cryogenic refrigerators used in existing plants are notdesigned to provide adequate cooling under such a pulsed regime.

Specifically, increasing the thermal load leads to an increase in theflow rate of the cooled helium, returned to the electromagnets of theplant.

This cools the entire refrigerator, because the cooling engendered bythe thermal load unbalances the exchangers 20, 21, 22, 23 and 24 betweenthe high-pressure section and low-pressure section. In addition, theevaporation caused by the load on the helium bath instantaneouslyincreases the return flow rate on the low-pressure side of each heatexchanger, thereby unbalancing the entire cryogenic refrigerator. Asubstantial increase in the thermal load may even cause the cryogenicrefrigerator to shut down.

To overcome this problem, it has already been suggested to smooth theimpact of the pulsed magnetic field on the variation of the thermal loadapplied to the cryogenic refrigerator.

This smoothing consists in limiting variations in the thermal loadactually applied to the cryogenic refrigerator, in order to ensure thenominal operation of the cryogenic refrigerator and therefore especiallyto prevent the refrigerator from shutting down.

For this purpose, it has been suggested to implement mechanical and/orthermal methods, either by installing dedicated means inside thecryogenic refrigerator, or by installing an additional device betweenthe plant employing strong magnetic fields and its cryogenicrefrigerator.

For example, the document by Dauguet et al., “Advances in CryogenicEngineering: Transactions of the Cryogenic Engineering Conference”, CECVol. 53, edited by J. G. Weisend II, pp. 564-569 proposes keeping thethermal load on the cryogenic refrigerator constant.

To do this, the cryogenic refrigerator comprises many heat exchangersand cryogenic valves, these valves being activated in order to keep thecooling power delivered by the refrigerator stable at a valuecorresponding to the average thermal load applied by the tokamak. Thus,the cryogenic refrigerator operates in a permanent or almost permanentregime, even when the tokamak is operating in a pulsed regime.

A similar solution is proposed in document WO 2009/024705.

According to another example, document “Design of the ITER-FEATcryoplant to achieve stable operation over a wide range of experimentalparameters and operation scenarios”, by Claudet et al., FusionEngineering and Design, 58-59 (2001), pp. 205-209, suggests installingan intermediate device between the tokamak and the cryogenicrefrigerator.

This device allows some of the additional helium flow obtained during apeak in the thermal load applied to the cryogenic refrigerator to bediverted upstream of the refrigerator. Thus, the cryogenic refrigeratordoes not see the increased helium flow related to the pulsed operatingregime of the tokamak.

These proposed solutions are based on mechanical and/or thermal devicesintended to smooth the variation in the thermal load liable to beapplied to the cryogenic refrigerator.

These solutions work correctly.

However, they require additional parts (heat exchangers, valves, etc.)which may quickly prove to be expensive and they are, sometimes,difficult to implement.

One objective of the invention is to solve at least one of the drawbacksof existing cryogenic refrigerators.

Another objective is to solve at least one of the drawbacks of existingcryogenic refrigerators when the plant to be cooled operates in avariable regime.

As was mentioned above, a variable regime may be encountered in manyfields. Therefore, use of the invention is not limited to plantsgenerating strong pulsed magnetic fields, such as a tokamak, but extendsto any plant requiring cryogenic refrigeration.

To achieve at least one of these objectives, the invention provides amethod for estimating a thermal load applied to a cryogenic refrigeratorcomprising a phase separator containing a bath of liquid refrigerant, inwhich this thermal load is estimated using a computer program, saidprogram being based on a mass balance performed on the phase separatorallowing the variation in the time derivative of the height of the bathof liquid refrigerant in the phase separator to be expressed.

The method according to the invention will possibly have other technicalfeatures, whether alone or in combination:

-   -   the thermal load applied to the cryogenic refrigerator is a        variable thermal load;    -   the variable thermal load applied to the cryogenic refrigerator        is pulsed;    -   a step is provided for determining appropriate variables of the        cryogenic refrigerator, representing the variation in the time        derivative of the height of the bath of liquid refrigerant in        the phase separator, before the step of estimating the thermal        load using the computer program;    -   the cryogenic refrigerator comprising a valve at the inlet of        the phase separator, said appropriate variables comprise at        least the degree of opening of the valve, the temperature        upstream of the valve and their respective time derivatives, and        the thermal load; and    -   the computer program comprises a step in which the data relating        to the variation in the time derivative of the height of the        bath of liquid refrigerant in the phase separator are filtered.

To achieve at least one of these objectives, the invention also providesa method for regulating the cryogenic refrigerator subjected to athermal load, in which the thermal load applied to the refrigerator isestimated using the method for estimating a thermal load applied to thisrefrigerator, according to the invention, and then at least oneoperating parameter of the refrigerator is regulated depending on thevalue of the thermal load estimated beforehand.

The regulating method according to the invention will possibly haveother technical features, in particular:

-   -   a step in which the regulation is implemented by modifying the        degree of opening of the valve of the refrigerator, said valve        being located at the inlet of the phase separator.

To achieve at least one of these objectives, the invention also providesa computer program product comprising programming code instructions forimplementing the method for estimating a thermal load applied to acryogenic refrigerator, according to the invention.

Other features, aims and advantages of the invention will becomeapparent from the following detailed description, given with referenceto the following figures:

FIG. 3 is a schematic showing a mass balance performed on the phaseseparator of the cryogenic refrigerator, for example as illustrated inFIG. 1 or 4;

FIG. 4 is a schematic of the cryogenic refrigerator named “station800W@4.5K” and installed at the Commissariat à l'énergie atomique, onwhich the method according to the invention was tested and validated;

FIG. 5 is a graph showing, on the one hand, the variation in the thermalload actually applied to the cryogenic refrigerator shown in FIG. 4during a trial, and on the other hand, the variation of the thermal loadestimated using the method according to the invention.

The method especially comprises a step of estimating the thermal load wapplied to the cryogenic refrigerator 1, 1′ using a computer program.

Details of the development of this program are given below.

This program especially operates on a mass conservation balanceperformed on the phase separator 40, 40′ of the cryogenic refrigerator1, 1′.

This mass balance performed on the phase separator 40, 40′ is shownschematically in FIG. 3. Typically, it is a phase separator 40, 40′ suchas shown in FIGS. 1 and 4, respectively, which comprises a bath 41, 41′of liquid refrigerant, for example of helium.

The mass conservation equation for the phase separator 40, 40′ iswritten:

$\begin{matrix}{\frac{h}{t} = {{m_{+} - m_{-}} = {{f\left( {u_{1},P,T} \right)} - \frac{w}{L_{v}}}}} & \left( {{Eq}.\mspace{14mu} 1} \right)\end{matrix}$

where:

h is the height of the bath 41, 41′ of liquid refrigerant in the phaseseparator 40, 40′ (% of the maximum attainable height h_(max) of thebath);

m₊ is the flow rate of gas entering into the phase separator (g/s);

m⁻ is the flow rate of gas leaving the phase separator (g/s);

w is the thermal load applied to the phase separator, i.e. to therefrigerator, by the plant (W);

u₁ is the degree of opening of the valve 50, 50′ (%);

L_(V) is the latent vaporisation heat (J/g); and

f(u₁, P, T) is a function depending on the degree of opening u₁ of thevalue 50, 50′, on the pressure P upstream of the valve and on thetemperature T upstream of this value.

The function f(u₁, P, T) may take various developed forms depending onthe precision desired for the estimation of the time derivative

$\frac{h}{t}$

of the height h of the liquid refrigerant in the phase separator 40,40′, which is denoted {dot over (h)} in the following.

In the example given below, the quantity {dot over (h)} depends on thefollowing parameters: u₁, T·u₁, {dot over (u)}₁, {dot over (T)} and w,where {dot over (u)}₁ and {dot over (T)} are the time derivatives of thedegree of opening of the valve 50, 50′ and of the temperature upstreamof this valve, respectively. These parameters represent the mainparameters influencing the mass balance performed on the phase separator40, 40′.

In this example, the influence of the pressure P in the function f isnot taken into consideration.

Different parameters may be chosen, depending on the desired precisionof the model, in particular parameters relating to the pressure P couldbe incorporated.

According to another example, the model takes into account, in order todefine the quantity {dot over (h)}, the following parameters: u₁, T, andtheir respective time derivatives, and the thermal load w.

It is then sought to express the equation (Eq. 1) in a linear form, i.e.in the form of (Eq. 2):

$\overset{.}{h} = {\left( {u_{1},\left( {Tu}_{1} \right),{\overset{.}{u}}_{1},\overset{.}{T},w} \right) \cdot \begin{pmatrix}a_{1} \\a_{2} \\a_{3} \\a_{4} \\a_{5}\end{pmatrix}}$

where: a₁, a₂, a₃, a₄ and a₅ are coefficients to be determined.

The equation (Eq. 2) may be written in the form of the equation (Eq. 3):

{dot over (h)}=L _(q) q+L _(z) ż+a ₅ w

on account of the following notations:

${q:=\begin{pmatrix}u_{1} \\{Tu}_{1}\end{pmatrix}};$ ${z = \begin{pmatrix}u_{1} \\T\end{pmatrix}};$ L_(q)(a₁, a₂); L_(z)(a₃, a₄)

To estimate the thermal load w applied to the cryogenic refrigerator 1,1′ by the plant, the quantity of interest is the filtered timederivative {dot over (h)} of the height h of liquid refrigerant in thephase separator, which is expressed in a Laplace field.

Filtering allows any “noise” that is capable of influencing thedetermination of the variable {dot over (h)} to be removed.

Thus, if the time derivative {dot over (h)} of this height h is denotedH(s) in the Laplace field, and the function obtained after first-orderfiltering of the quantity H(s) is denoted Y(s), then a priori theirrelationship can be expressed in the form of the equation (Eq. 4):

${Y(s)} = {\frac{s}{1 + {\tau_{f}s}} \cdot {H(s)}}$

where τ_(f) is the time constant of this first-order filter, which mustbe defined.

The time constant is defined in the following way. The valve 50, 50′ isopened with a given opening profile, a stepped profile for example.Then, the variation of the time derivative of the height of the bath ofliquid is monitored. It is then possible to determine the time constantin a way known per se to those skilled in the art.

Equation (Eq. 4) then corresponds, in the real field, to the followingdifferential equation (Eq. 5):

$\overset{.}{y} = {\frac{1}{\tau_{f}}\left\lbrack {{- y} + \overset{.}{h}} \right\rbrack}$

where: {dot over (y)} is the real variable associated with the filteredfunction Y(s) of equation (Eq. 4) in the Laplace field, and y is thereal integrated value of the variable {dot over (y)}.

By inserting the relationship given by equation (Eq. 3) obtained fromthe mass conservation balance performed on the phase separator 40, 40′into equation (Eq. 5) above, equation (Eq. 6) is obtained:

$\overset{.}{y} = {\frac{1}{\tau_{f}}\left\lbrack {{- y} + {L_{q}q} + {L_{z}\overset{.}{z}} + {a_{5}w}} \right\rbrack}$

It will be understood that this differential equation expresses the massconservation balance performed on the phase separator after first-orderfiltering has been carried out on the variable {dot over (h)}.

Next, the state vector x_(w) is introduced in the form of the equality(Eq. 7):

x _(w):=τ_(f) ·y−L _(z) ·z

This notation then allows equation (Eq. 6) to be written in the form ofthe following equation (Eq. 8):

{dot over (x)} _(w) =−y+L _(q) q+a ₅ w

Next, by replacing the variable y with the variable

$\frac{1}{\tau_{f}}\left\lbrack {x_{w} + {L_{z}z}} \right\rbrack$

the following equation is obtained:

${\overset{.}{x}}_{w} = {{- {\frac{1}{\tau_{f}}\left\lbrack {x_{w} + {L_{z}z}} \right\rbrack}} + {L_{q}q} + {a_{5}w}}$

which finally gives the following state representation (Eq. 9):

${\overset{.}{x}}_{w} = {{\left\lbrack \frac{- 1}{\tau_{f}} \right\rbrack x_{w}} + {\left( {L_{q},\frac{- L_{z}}{\tau_{f}}} \right)\begin{pmatrix}q \\z\end{pmatrix}} + {a_{5}w}}$$y = {{\left\lbrack \frac{1}{\tau_{f}} \right\rbrack x_{w}} + {\left\lbrack \left( {0,\frac{L_{z}}{\tau_{f}}} \right) \right\rbrack \begin{pmatrix}q \\z\end{pmatrix}}}$

where x_(w) is the state and q, z and w the inputs and y the output.

For the sake of clarity, the state representation (Eq. 9) may also beexpressed in the following form (Eq. 10):

${\overset{.}{x}}_{w} = {{Ax}_{w} + {B\begin{pmatrix}q \\z\end{pmatrix}} + {Gw}}$ $y = {{Cx}_{w} + {D\begin{pmatrix}q \\z\end{pmatrix}}}$

on account of the following matrix notations:

${A = \left\lbrack \frac{- 1}{\tau_{f}} \right\rbrack};$${B = \left( {L_{q},\frac{- L_{z}}{\tau_{f}}} \right)};$${C = \left\lbrack \frac{1}{\tau_{f}} \right\rbrack};$${D = \left\lbrack \left( {0,\frac{L_{z}}{\tau_{f}}} \right) \right\rbrack};$G = a₅

The state representation (Eq. 9/Eq. 10) finally allows the filteredvariable y representative of the height h of the liquid refrigerant inthe phase separator to be expressed as a function of q, z, w and thecoefficients a₁, a₂, a₃, a₄, a₅ and τ_(f).

This state representation can be easily implemented in a programmablecontroller.

This state representation may be extended to estimate the thermal loadapplied to the cryogenic refrigerator.

To do this, the extended state of the state representation (Eq. 9) isdefined by the following matrix, (Eq. 11):

$\xi:=\begin{pmatrix}x_{w} \\s_{w}\end{pmatrix}$

where the second state vector s_(w) represents the thermal load w.

Then, the state representation (Eq. 10) can be expressed in thefollowing form, denoted (Eq. 12):

$\overset{.}{\xi}:={\begin{pmatrix}{\overset{.}{x}}_{w} \\{\overset{.}{s}}_{w}\end{pmatrix} = {{\underset{\underset{\overset{\_}{A}}{}}{\begin{pmatrix}A & G \\0 & 0\end{pmatrix}}\begin{pmatrix}x_{w} \\s_{w}\end{pmatrix}} + {\begin{pmatrix}B \\0_{1 \times 4}\end{pmatrix} \cdot \begin{pmatrix}q \\z\end{pmatrix}}}}$$y = {{\underset{\underset{\overset{\_}{C}}{}}{\begin{pmatrix}C & 0\end{pmatrix}}\xi} + {D\begin{pmatrix}q \\z\end{pmatrix}}}$

This representation assumes that {dot over (s)}_(w)=0, i.e. that thevariable s_(w) representative of the thermal load w applied to therefrigerator has a constant profile. This assumption is correct when theperturbation that arrives at the bath of liquid takes a step form. Thisis especially the case for a tokamak.

Of course, other assumptions could be made. For example, it may beassumed that the perturbation arrives in the form of a ramp or asinusoidal variation. This could be the case for cryogenic refrigeratorsused in plants other than tokamaks.

If (Eq. 13) is written:

$\begin{pmatrix}q \\z\end{pmatrix}:={\begin{pmatrix}u_{1} \\{Tu}_{1} \\u_{1} \\T\end{pmatrix}\mspace{14mu} \text{=:}\mspace{11mu} {U\left( {T,u_{1}} \right)}}$

then the state representation of (Eq. 12) can be written in the form of(Eq. 14):

{dot over (ξ)}=[Ā]·ξ+[ B]·U(Y,u ₁)

y=[ C]·ξ+[D]·U(T,u ₁)

Next, if discrete matrices representative of the matrices Ā and B for asampling period T_(e) are denoted A_(d) and B_(d) respectively, theabove relationships can be written:

ξ⁺ =[A _(d) ]·ξ+[B _(d) ]·U(T,u ₁)

y=[C]·ξ+[D]·U(T,u ₁)

where ξ⁺ is the discrete form of the matrix ξ.

It will be recalled that when the continuous matrices A, B, C, D, of astate representation are known, the associated discrete matrices A_(d),B_(d), C_(d) and D_(d) are written:

-   A_(d)=e^(A·T); where T_(e) is the sampling period;

B_(d) = ∫₀^(T_(e))^(A ⋅ s) ⋅ B ⋅ s;

-   C_(d)=C; and-   D_(d)=D.

The relationships expressing ξ⁺ and y can be written in shortened form.

To do this, the observer equation will be recalled, namely:

{circumflex over (ξ)}⁺ =A _(d) {circumflex over (ξ)}+B _(d) U+L(y−ŷ)

where L is a two-row column matrix comprising real values called theobserver gain. The observer gain only depends on the response time t_(r)chosen for the observer.

Next, the relationship expressing y is used for the variable ŷ. In thiscase ŷ can be written:

ŷ=[ C]·{circumflex over (ξ)}+[D]·U(T,u ₁).

By replacing ŷ with the latter expression in the observer equation, theshortened form mentioned above is finally obtained, namely equation (Eq.15):

{circumflex over (ξ)}⁺ =[A _(d) −L C]·{circumflex over (ξ)}+[B _(d)−LD]·U(T,u ₁)+L·y

in which the observer gain L is chosen so that the eigenvalues of thematrix A_(d)−L C are in part real and strictly negative. As should beunderstood with ease by those skilled in the art, the eigenvalues of thematrix A_(d)−L C correspond to the poles of the corresponding transferfunction.

These poles are especially related to the time constant τ_(f), insofaras this parameter appears in the expression of this matrix. Thus, thesepoles are related to the response time t_(r) chosen for the observer.

However, it must be noted that the variable y is not present in equation(Eq. 15). It must therefore be determined in some other way. To do this,the differential equation (Eq. 5) is used, namely:

$\overset{.}{y} = {{\frac{1}{\tau_{f}}\left\lbrack {{- y} + \overset{.}{h}} \right\rbrack}.}$

Next, the vector:

η:=τ_(f) ·y−h

is defined and equation (Eq. 5) is rewritten in the following form:

$\overset{.}{\eta} = {{- y} = {\frac{- 1}{\tau_{f}}\left\lbrack {\eta + h} \right\rbrack}}$$y = {\frac{1}{\tau_{f}}\left\lbrack {\eta + h} \right\rbrack}$

which may be expressed in the form of the following staterepresentation:

$\overset{.}{\eta} = {{\left\lbrack \frac{- 1}{\tau_{f}} \right\rbrack \cdot \eta} + {\left\lbrack \frac{- 1}{\tau_{f}} \right\rbrack \cdot h}}$$y = {{\left\lbrack \frac{1}{\tau_{f}} \right\rbrack \cdot \eta} + {\left\lbrack \frac{1}{\tau_{f}} \right\rbrack \cdot {h.}}}$

This state representation may then be expressed in discrete form, for asampling period T_(e), in the following way (Eq. 16):

η⁺ = [R_(d)] ⋅ η + [E_(d)] ⋅ h$y = {{\left\lbrack \frac{1}{\tau_{f}} \right\rbrack \cdot \eta} + {\left\lbrack \frac{1}{\tau_{f}} \right\rbrack \cdot h}}$

where η⁺ is the discrete form of the variable {dot over (η)}, thematrices [R_(d)] and [E_(d)] being the discrete matrices of the staterepresentation of equation (Eq. 5). These matrices are defined, as afunction of the continuous matrices [R]=[−1/τ_(f)] and E=[−1/τ_(f)], bythe following relationships:

-   R_(d)=e^(R·T) ^(e) where T_(e) is the sampling period; and

E_(d) = ∫₀^(T_(e))^(R ⋅ s) ⋅ E ⋅ s.

By then grouping equations (Eq. 15) and (Eq. 16), a set of equations(Eq. 17) may be written:

     η⁺ = [R_(d)] ⋅ η + [E_(d)] ⋅ h${\hat{\xi}}^{+} = {{{\left\lbrack {A_{d} - {L\overset{\_}{C}}} \right\rbrack \cdot \hat{\xi}} + {\left\lbrack {B_{d} - {LD}} \right\rbrack \cdot {U\left( {T,u_{1}} \right)}} + {L \cdot {\left( {{\left\lbrack \frac{1}{\tau_{f}} \right\rbrack \cdot \eta} + {\left\lbrack \frac{1}{\tau_{f}} \right\rbrack \cdot h}} \right).\mspace{79mu} \hat{w}}}} = {\left( {0,1} \right)\hat{\xi}}}$

This set of equations may be expressed in the following form (Eq. 18):

$X^{+} = {{\begin{pmatrix}R_{d} & 0_{1 \times 2} \\{\frac{1}{\tau_{f}}L} & {A_{d} - {L\overset{\_}{C}}}\end{pmatrix}X} + {\begin{pmatrix}0_{1 \times 4} & E_{d} \\{B_{d} - {LD}} & {\frac{1}{\tau_{f}}L}\end{pmatrix}\begin{pmatrix}{U\left( {T,u_{1}} \right)} \\h\end{pmatrix}}}$ ŵ = (0, 0, 1)X

if the following notation

$X:=\begin{pmatrix}\eta \\\hat{\xi}\end{pmatrix}$

is used; or, more simply (Eq. 19):

X ⁺ =A _(obs) X+B _(obs) U _(obs)(T,u ₁ ,h)

ŵ=C_(obs)X

with the following matrices:

${A_{obs} = \begin{pmatrix}R_{d} & 0_{1 \times 2} \\{\frac{1}{\tau_{f}}L} & {A_{d} - {L\overset{\_}{C}}}\end{pmatrix}};$ ${B_{obs} = \begin{pmatrix}0_{1 \times 4} & E_{d} \\{B_{d} - {LD}} & {\frac{1}{\tau_{f}}L}\end{pmatrix}};$ C_(obs) = (0, 0, 1); and${U_{obs}\left( {T,u_{1},h} \right)} = {\begin{pmatrix}u_{1} \\{Tu}_{1} \\u_{1} \\T \\h\end{pmatrix}.}$

Equation (Eq. 19) thus allows the thermal load w applied to therefrigerator to be estimated, in a form that can be easily implementedin a programmable controller.

However the thermal load w can be estimated, using (Eq. 19), only oncethe matrices A_(obs) and B_(obs) have been completely defined, whichrequires that the coefficients a₁ to a₅ be identified, the values of thematrix L be determined and the time constant τ_(f) of the first-orderfilter be chosen.

The coefficients a₁ to a₅ are identified in the following way.

The equation (Eq. 2) is integrated between a reference time t₁ and thetime t, and then filtered.

The integration returns equation (Eq. 20):

h(t)−h(t ₁)=M(t)*a

where M(t) is the line matrix given by:

${{M(t)} = \begin{pmatrix}{\int_{t_{1}}^{t}{{u(\tau)}{\tau}}} & {\int_{t_{1}}^{t}{{T(\tau)}{u_{1}(\tau)}{\tau}}} & {{u_{1}(t)} - {u\left( t_{1} \right)}} & {{T(t)} - {T\left( t_{1} \right)}} & {\int_{t_{1}}^{t}{{w(\tau)}{\tau}}}\end{pmatrix}},$

and a is the column vector of the coefficients a₁ to a₅.

Equation (Eq. 20) is then filtered with the following filter:

${F(s)} = {\frac{s}{1 + {\tau_{f}s}}.}$

The coefficients a₁ to a₅ are then determined by minimizing, using theleast squares method, the relationships forming a system of equations(Eq. 21), the unknowns of which are the coefficients of the matrix a:

${{\begin{pmatrix}{{F(M)}\left( t_{1} \right)} \\{{F(M)}\left( t_{2} \right)} \\\vdots \\{{F(M)}\left( t_{N} \right)}\end{pmatrix} \cdot a} - \begin{pmatrix}0 \\{{{F(h)}\left( t_{2} \right)} - {h\left( t_{1} \right)}} \\\vdots \\{{{F(h)}\left( t_{N} \right)} - {h\left( t_{1} \right)}}\end{pmatrix}}$

which allows equation (Eq. 20) to be solved in its filtered form betweenthe initial time t₁ and the final time t_(N), where N=5 in thisinstance.

To do this, it is necessary to measure the values of the coefficients ofthe other matrices.

These measurements were carried out on the refrigerator 1′ shownschematically in FIG. 4.

This refrigerator 1′ is similar to the refrigerator 1 described withreference to FIG. 1, especially as regards the phase separator 40′. Thisrefrigerator 1′ however differs slightly from the refrigerator describedwith reference to FIG. 1 in various ways.

Specifically, it only has four heat exchangers 20′, 21′, 22′, 23′ forlowering the temperature of the refrigerant, in this case helium, from300 K to 4.5 K. The pressure delivered by the compressor is 16 bar.Moreover, only one turbine 31′ is provided and the first heat exchanger20′ contains an additional heat exchanger 60′ (nitrogen, 80 K).

The coefficients of the matrix a were calculated for a compressor outletpressure P of 16 bar. Moreover, the temperature upstream of this valve50′ is T=7 K. The time t₁ considered is t₁=1 s, at which point datarecording begins. The interval between measurements is 3 s, i.e.dt=t_(i+1)−t_(i)=3 s, for i ranging from 1 to N.

Under these conditions, solving the system of equations (Eq. 21) returnsthe coefficients in table 1 below:

TABLE 1 Coefficient of the matrix a value a₁ 0.003710768747845 a₂−0.000396506514909 a₃ 0.024617305832635 a₄ 0.286747643082885 a₅−0.000240828045189

Once the coefficients a₁ to a₅ have been obtained, they are entered intothe program, as are the time constant τ_(f) and the coefficients of theobserver gain L.

A 5% response time, denoted tr5%, was chosen for the observation. Inthis instance, tr5%=60 s. This indicates that the observer is delayed inorder to obtain a 5% estimate of the thermal load w in 60 s.

The time constant τ_(f) is determined using the relationshiptr5%=3τ_(f). Thus, the time constant τ_(f) has a value τ_(f)=20 s.

Moreover, to determine the eigenvalues of the matrix A_(d)−L C, thepoles z of the corresponding transfer function are determined, i.e.

$z = {^{- \frac{3 \cdot {t}}{{tr}\; 5\%}} = \begin{pmatrix}0.8607 \\0.8607\end{pmatrix}}$

insofar as dt=3 s, as mentioned above.

These data especially allow the matrices A_(obs) and B_(obs) to bedefined.

The thermal load w can then be estimated with this program using thesystem of equations (Eq. 19).

FIG. 5 shows the estimation of the thermal load obtained by the programdescribed above and the actual variation in the thermal load applied tothe refrigerator shown in FIG. 4.

As may be seen in said FIG. 5, although the thermal load applied to therefrigerator 1′ is highly variable, and even random, the estimation isof a very high quality. The estimation follows the actual variation inthe thermal load during the peaks in power.

This estimation is also very good in an equilibrium regime, for examplebetween the time t=0 s and the time t=3×10⁴ s in FIG. 5.

Finally, it is possible to correctly estimate, in real time, thevariation in the thermal load applied to the cryogenic refrigerator.

This estimation is very good whether the operating regime is a variableor equilibrium regime.

This is particularly advantageous because it is then possible toregulate one or more operating parameters of the cryogenic refrigerator.

It may for example be envisioned to regulate the temperature at theoutlet of the turbine 31′, the height of the bath of liquid in theseparator, or other parameters.

Considering the above example, it will then be possible to physicallyadjust variables such as the degree of opening of the valve, thetemperature and/or the pressure upstream of this valve, etc. in order toregulate the or each operating parameter considered.

In particular, in the case of pulsed tokamak regimes, it is possible toestimate, in real time, the value of the thermal load applied to thecryogenic refrigerator, without knowledge of events outside of therefrigerator.

This regulation then allows adequate cooling of the plant to be ensured,avoiding the risk of the refrigerator shutting down and consequently theplant itself shutting down.

This regulation is inexpensive because it requires no major hardware.

It should be noted that the coefficients identified for the matrix a arevalid for the refrigerator 1′ shown in FIG. 4, taking account of theparameters chosen for equation (Eq. 2).

Different values would have to be identified for the matrix a if otherparameters were to be included in equations (Eq. 2). For example, if itwere desired to have an even more precise model, it could be envisionedto furthermore take into account the pressure P upstream of the valve50′, and the time derivative of this pressure P.

Different values would also have to be identified for the matrix a forthe same refrigerator operating under different conditions, for exampleif the outlet pressure of the compressor were different.

1. A method for estimating a thermal load applied to a cryogenicrefrigerator comprising a phase separator containing a bath of liquidrefrigerant, in which this thermal load is estimated using a computerprogram, said program being based on a mass balance performed on thephase separator allowing the variation in the time derivative of theheight of the bath of liquid refrigerant in the phase separator to beexpressed.
 2. The method as claimed in claim 1, in which the thermalload applied to the cryogenic refrigerator is a variable thermal load.3. The method as claimed in claim 1, in which the variable thermal loadapplied to the cryogenic refrigerator is pulsed.
 4. The method asclaimed in claim 1, in which a step is provided for determiningappropriate variables of the cryogenic refrigerator, representing thevariation in the time derivative of the height of the bath of liquidrefrigerant in the phase separator, before the step of estimating thethermal load using the computer program.
 5. The method as claimed inclaim 1, in which, the cryogenic refrigerator comprising a valve at theinlet of the phase separator, said appropriate variables comprise atleast the degree of opening of the valve, the temperature upstream ofthe valve and their respective time derivatives, and the thermal load.6. The method as claimed in claim 4, in which the computer programcomprises a step in which the data relating to the variation in the timederivative of the height of the bath of liquid refrigerant in the phaseseparator are filtered.
 7. A method for regulating the cryogenicrefrigerator subjected to a thermal load, in which the thermal loadapplied to the refrigerator is estimated using the method as claimed inclaim 1, and then at least one operating parameter of the refrigeratoris regulated depending on the value of the thermal load estimatedbeforehand.
 8. The method as claimed in claim 7, in which the regulatingstep is implemented by modifying the degree of opening of the valve ofthe refrigerator, said valve being located at the inlet of the phaseseparator.
 9. A computer program product comprising programming codeinstructions for implementing a method as claimed in claim 1.